home *** CD-ROM | disk | FTP | other *** search
/ Skunkware 5 / Skunkware 5.iso / src / X11 / endo / maps.c < prev    next >
C/C++ Source or Header  |  1995-05-03  |  22KB  |  1,083 lines
  1. /*************************************************************************
  2.  *                                                                       *
  3.  *  Copyright (c) 1992, 1993 Ronald Joe Record                           *
  4.  *                                                                       *
  5.  *  All rights reserved. No part of this program or publication may be   *
  6.  *  reproduced, transmitted, transcribed, stored in a retrieval system,  *
  7.  *  or translated into any language or computer language, in any form or *
  8.  *  by any means, electronic, mechanical, magnetic, optical, chemical,   *
  9.  *  biological, or otherwise, without the prior written permission of:   *
  10.  *                                                                       *
  11.  *      Ronald Joe Record (408) 458-3718                                 *
  12.  *      212 Owen St., Santa Cruz, California 95062 USA                   *
  13.  *                                                                       *
  14.  *************************************************************************/
  15.  
  16. #include <math.h>
  17. #ifndef NeXT
  18. #include <values.h>
  19. #endif
  20. #include "defines.h"
  21.  
  22. double
  23. wrap(z, a, b)
  24. double z;
  25. double a, b;
  26. {
  27.     double w;
  28.  
  29.     w = z;
  30.     while (w > b) 
  31.         w -= b - a;
  32.     while (w < a) 
  33.         w += b - a;
  34.     return(w);
  35. }
  36.  
  37. double
  38. trace(x, y, p)
  39. double x, y;
  40. double *p;
  41. {
  42.     static pair c, cx, cy;
  43.     static double r;
  44.     extern PFP map;
  45.     extern double delta;
  46.  
  47.     c = (*map)(x, y, p);
  48.     cx = (*map)(x+delta, y, p);
  49.     cy = (*map)(x, y+delta, p);
  50.  
  51.     r = (cx.x + cy.y - c.x - c.y) / delta;
  52.     return(r);
  53. }
  54.  
  55. pair
  56. gaertner(x, y, p)        /* the Gaertner "logistic" map */
  57. double x, y; 
  58. double *p;
  59. {
  60.     pair coord;
  61.  
  62.     coord.x=(p[0]*x*(p[4] - (p[5]*x))) + (p[1]*y*(p[6] - (p[7]*y)));
  63.     coord.y=(p[2]*y*(p[8] - (p[9]*y))) + (p[3]*x*(p[10] - (p[11]*x)));
  64.     return(coord);
  65. }
  66.  
  67. pair
  68. dgaertner(x, y, p)
  69. double x, y;
  70. double *p;
  71. {
  72.     pair r;
  73.     static double AC, BD;
  74.  
  75.     AC = p[0] * p[2]; BD = p[1] * p[3];
  76.     /* calculate the determinant */
  77.     r.x = (AC*p[4]*p[8])-(BD*p[6]*p[10]);
  78.     r.x -= 2.0*((AC*p[5]*p[8])+(BD*p[6]*p[11]))*x;
  79.     r.x -= 2.0*((AC*p[4]*p[9]) + (BD*p[7]*p[10]))*y;
  80.     r.x += 4.0*((AC*p[5]*p[9]) - (BD*p[7]*p[11]))*x*y;
  81.     /* calculate the trace */
  82.     r.y = (p[0]*p[4]) + (p[2]*p[8]) - (2.0*((p[0]*p[5]*x) + (p[2]*p[9]*y)));
  83.     return(r);
  84. }
  85.  
  86. pair
  87. goodwin(x, y, p)            /* the Goodwin map */
  88. double x, y; 
  89. double *p;
  90. {
  91.     pair coord;
  92.  
  93.     if (p[0] == -x) {
  94.         coord.y = y;
  95.         coord.x = x;
  96.         return(coord);
  97.     }
  98.     coord.y = ((1.0 / (p[0] + x)) - (p[3] + p[4])) * y;
  99.     if (coord.y == 1.0) {
  100.         coord.y = y;
  101.         coord.x = x;
  102.         return(coord);
  103.     }
  104.     coord.x = (1.0 + p[1] * ((1.0 / (1.0 - coord.y)) - p[2]) - p[3]) * x;
  105.     return(coord);
  106. }
  107.  
  108. pair
  109. dgoodwin(x, y, p)
  110. double x, y;
  111. double *p;
  112. {
  113.     pair r;
  114.     double s;
  115.     static double A, B, C, D, E, yn1, sqy, aplusx, sqa;
  116.  
  117.     r.x = r.y = MAXDOUBLE;
  118.     A = p[0]; B = p[1]; C = p[2]; D = p[3]; E = p[4];
  119.     /* calculate the determinant */
  120.     if (A == -x)
  121.         return(r);
  122.     yn1 = 1.0 - (y/(A+x)) + (D*y) + (E*y); aplusx = A + x;
  123.     if ((aplusx == 0.0) || (yn1 == 0.0))
  124.         return(r);
  125.     sqy = yn1 * yn1; sqa = aplusx * aplusx;
  126.     r.x = 1.0 + ((((B*yn1) - (B*x*y/sqa)))/sqy) -(B*C) - D;
  127.     r.x = r.x * ((1.0/aplusx) - D - E);
  128.     s = (y/sqa)*(B*x*((-1.0/aplusx) + D + E))/sqy;
  129.     r.x = r.x - s;
  130.     /* calculate the trace later */
  131.     r.y = trace(x, y, p);
  132.     return(r);
  133. }
  134.  
  135. pair
  136. gucken(x, y, p)        /* the Guckenhiemer, Oster, Ipaktchi map */
  137. double x, y; 
  138. double *p;
  139. {
  140.     pair coord;
  141.     extern double exp();
  142.  
  143.     coord.x = ((p[0] * x) + (p[1] * y)) * exp(-p[2] * (x + y));
  144.     coord.y = p[3] * x;
  145.     return(coord);
  146. }
  147.  
  148. pair
  149. dgucken(x, y, p)
  150. double x, y;
  151. double *p;
  152. {
  153.     static pair r;
  154.     static double s, t;
  155.     extern double exp();
  156.  
  157.     s = p[2]*((p[0]*x) + (p[1]*y));
  158.     t = exp(-p[2]*(x+y));
  159.     /* calculate the determinant */
  160.     r.x = -p[3] * ((p[1] - s) * t);
  161.     /* calculate the trace */
  162.     r.y = (p[0] - s) * t;
  163.     return(r);
  164. }
  165.  
  166. pair
  167. standard(x, y, p) /* the Standard map (see Wood, Lichtenberg & Lieberman) */
  168. double x, y; 
  169. double *p;
  170. {
  171.     static pair coord;
  172.     static double s;
  173.     extern double sin();
  174.  
  175.     s = sin(y);
  176.     coord.x = wrap(x + (p[0] * s), 0.0, M_2PI);
  177.     coord.y = wrap(y + x + (p[1] * s), 0.0, M_2PI);
  178.     return(coord);
  179. }
  180.  
  181. pair
  182. dstandard(x, y, p)
  183. double x, y;
  184. double *p;
  185. {
  186.     pair r;
  187.     extern double cos();
  188.  
  189.     /* calculate the determinant */
  190.     r.x = 1.0 + ((p[1] - p[0]) * cos(y));
  191.     /* calculate the trace later */
  192.     r.y = trace(x, y, p);
  193.     return(r);
  194. }
  195.  
  196. pair
  197. hump(x, y, p)        /* the "double sine hump" map on the plane */
  198. double x, y; 
  199. double *p;
  200. {
  201.     pair coord;
  202.     static double sx, sy;
  203.     extern double sin();
  204.  
  205.     sx = sin(M_PI * x);
  206.     sy = sin(M_PI * y);
  207.     coord.x = wrap((p[0] * sx) + (p[1] * sy), 0.0, 2.0);
  208.     coord.y = wrap((p[2] * sy) + (p[3] * sx), 0.0, 2.0);
  209.     return(coord);
  210. }
  211.  
  212. pair
  213. dhump(x, y, p)
  214. double x, y;
  215. double *p;
  216. {
  217.     pair r;
  218.     extern double cos();
  219.  
  220.     /* calculate the determinant */
  221.     r.x = ((p[0]*p[2])-(p[1]*p[3])) * (M_PI*M_PI*cos(M_PI*x)*cos(M_PI*y));
  222.     /* calculate the trace later */
  223.     r.y = trace(x, y, p);
  224.     return(r);
  225. }
  226.  
  227. pair
  228. circle(x, y, p)        /* the "double circle" map on the plane */
  229. double x, y; 
  230. double *p;
  231. {
  232.     static pair c;
  233.     static double sx, sy;
  234.     extern double sin();
  235.  
  236.     sx = sin(M_2PI * x);
  237.     sy = sin(M_2PI * y);
  238.     c.x=wrap((p[4]+(p[8]*x)-(p[0]*sx))+(p[5]+(p[9]*y)-(p[1]*sy)),0.0,1.0);
  239.     c.y=wrap((p[6]+(p[10]*y)-(p[2]*sy))+(p[7]+(p[11]*x)-(p[3]*sx)),0.0,1.0);
  240.     return(c);
  241. }
  242.  
  243. pair
  244. dcircle(x, y, p)
  245. double x, y;
  246. double *p;
  247. {
  248.     pair r;
  249.     static double cy, cx;
  250.     extern double cos();
  251.  
  252.     cy = cos(M_2PI*wrap(y, 0.0, 1.0)); cx = cos(M_2PI*wrap(x, 0.0, 1.0));
  253.     /* calculate the determinant */
  254.     r.x = ((p[8]*p[10])-(p[9]*p[11]))+(((p[1]*p[11])-(p[2]*p[8]))*M_2PI*cy);
  255.     r.x += ((p[3]*p[9]) - (p[0]*p[10])) * M_2PI * cx;
  256.     r.x += ((p[0]*p[2])-(p[1]*p[3]))*4.0*M_2PI*M_2PI* cx * cy;
  257.     /* calculate the trace later */
  258.     r.y = trace(x, y, p);
  259.     return(r);
  260. }
  261.  
  262. pair
  263. circle2(x, y, p)    /* a 2nd "double circle" map on the plane */
  264. double x, y; 
  265. double *p;
  266. {
  267.     pair coord;
  268.     extern double sin(), cos();
  269.  
  270.     coord.x = wrap((p[4]+(p[8]*x)-(p[0]*sin(M_2PI*x))) +
  271.             (p[5]+(p[9]*y)-(p[1]*cos(M_2PI*y))), 0.0, 1.0);
  272.     coord.y = wrap((p[6]+(p[10]*y)-(p[2]*sin(M_2PI*y))) +
  273.             (p[7]+(p[11]*x)-(p[3]*cos(M_2PI*x))), 0.0, 1.0);
  274.     return(coord);
  275. }
  276.  
  277. pair
  278. dcircle2(x, y, p)
  279. double x, y;
  280. double *p;
  281. {
  282.     pair r;
  283.     static double cy, cx, sy, sx;
  284.     extern double cos(), sin();
  285.  
  286.     cy = cos(M_2PI*y); cx = cos(M_2PI*x);
  287.     sy = sin(M_2PI*y); sx = sin(M_2PI*x);
  288.     /* calculate the determinant */
  289.     r.x = (p[9]*p[11])-(p[8]*p[10])+(p[1]*p[11]*M_2PI*sy)+(p[2]*p[8]*M_2PI*cy);
  290.     r.x += (p[3]*p[9]*M_2PI*sx) + (p[0]*p[10]*M_2PI*cx);
  291.     r.x += (p[1]*p[3]*4.0*M_2PI*M_2PI*sx*sy)-(p[0]*p[2]*4.0*M_2PI*M_2PI*cx*cy);
  292.     /* calculate the trace later */
  293.     r.y = trace(x, y, p);
  294.     return(r);
  295. }
  296.  
  297. #ifdef Gardini
  298.  
  299. pair
  300. gard(x, y, p) /* Gardini figure 8 from Barugola, A. & Cathala, J. C. */
  301. double x, y;   /* "Annular Chaotic Areas", Nonlinear Analysis, Theory, */
  302. double *p;     /* Methods, & Applications Vol 10, no 11, pp. 1223-1236, 1986 */
  303. {
  304.     static pair coord;
  305.  
  306.     coord.x = y - (p[1]*x*y);    /* p[1] = 2 */
  307.     coord.y = (y - (p[0]*x)) + (p[2]*x*x*x); /* p[2] = 2 */
  308.     return(coord);
  309. }
  310.  
  311. pair
  312. dgard(x, y, p)
  313. double x, y;
  314. double *p;
  315. {
  316.     static pair r;
  317.  
  318.     r.x = (3.0*p[1]*p[2]*x*x*x) - (p[1]*p[0]*x) - (p[1]*y);
  319.     r.x += (3.0*p[2]*x*x) -p[0];
  320.     r.y = 1.0 - (p[1]*y);
  321.     return(r);
  322. }
  323.  
  324. pair
  325. gard2(x, y, p) /* Gardini fig 18 from Cathala, J.C. "On Some Properties */
  326. double x, y;    /* of Absorbtive Areas in Second Order Endomorphisms" in the */
  327. double *p;      /* 1989 European Conference on Iteration Theory, pp. 42-54 */
  328. {
  329.     static pair coord;
  330.  
  331.     coord.x = (p[1]*x) + y;        /* 0 < b < 0.01 */
  332.     coord.y = p[0] + (x*x);        /* a = -0.7 */
  333.     return(coord);
  334. }
  335.  
  336. pair
  337. dgard2(x, y, p)
  338. double x, y;
  339. double *p;
  340. {
  341.     static pair r;
  342.  
  343.     r.x = -2.0 * x;
  344.     r.y = p[1];
  345.     return(r);
  346. }
  347.  
  348. pair
  349. gard3(x, y, p) /* Gardini figure 20 unknown source */
  350. double x, y; 
  351. double *p;
  352. {
  353.     static pair coord;
  354.  
  355.     coord.x = p[1] * y;    /* p[1] = 1 */
  356.     coord.y = y - (p[0]*x) + (x*x);
  357.     return(coord);
  358. }
  359.  
  360. pair
  361. dgard3(x, y, p)
  362. double x, y;
  363. double *p;
  364. {
  365.     static pair r;
  366.     
  367.     r.x = (p[0] - (2.0*x))*p[1];
  368.     r.y = 1.0;
  369.     return(r);
  370. }
  371.  
  372. #endif /* Gardini */
  373.  
  374. double
  375. poly(x, p)
  376. double x;
  377. double *p;
  378. {
  379.     static int i;
  380.     static double r, s;
  381.  
  382.     s = 1.0;
  383.     r = 0.0;
  384.     for  (i=0; i<p[0]; i++) {
  385.         r += p[i+1]*s;
  386.         s *= x;
  387.     }
  388.     return(r);
  389. }
  390.  
  391. double
  392. dpoly(x, p)
  393. double x;
  394. double *p;
  395. {
  396.     static int i;
  397.     static double r, s;
  398.  
  399.     s = 1.0;
  400.     r = 0.0;
  401.     for  (i=1; i<p[0]; i++) {
  402.         r += i*p[i+1]*s;
  403.         s *= x;
  404.     }
  405.     return(r);
  406. }
  407.  
  408. pair
  409. alexander(x, y, p)    /* Alexander map described in Science News 11/14/92 */
  410. double x, y;
  411. double *p;
  412. {
  413.     pair z;
  414.  
  415.     z.x = (x*x) - (y*y) - x - (p[0]*y);
  416.     z.y = (2.0*x*y) - (p[1]*x) + y;
  417.     return(z);
  418. }
  419.  
  420. pair
  421. dalexander(x, y, p)
  422. double x, y;
  423. double *p;
  424. {
  425.     pair r;
  426.  
  427.     r.x = (4.0*((x*x) + (y*y))) + ((p[0]-p[1])*2.0*y) - (p[0]*p[1]) - 1.0;
  428.     r.y = (4.0*x);
  429.     return(r);
  430. }
  431.  
  432. pair
  433. secant(x, y, p)    /* secant method for f(x) = ax^n + bx^n-1 + ... + px + q */
  434. double x, y;
  435. double *p;
  436. {
  437.     pair z;
  438.     static double fx, fy;
  439.     extern double delta;
  440.  
  441.     fx = poly(x, p);
  442.     fy = poly(y, p);
  443.     if (ABS(fx) < delta)
  444.         z.x = x;
  445.     else {
  446.         if (fx != fy)
  447.             z.x = x - ((fx * (x - y)) / (fx - fy));
  448.         else
  449.             z.x = MAXDOUBLE;
  450.     }
  451.     z.y = x;
  452.     return(z);
  453. }
  454.  
  455. pair
  456. dsecant(x, y, p)
  457. double x, y;
  458. double *p;
  459. {
  460.     static double fx, fy, s;
  461.     static pair r;
  462.  
  463.     fx = poly(x, p);
  464.     fy = poly(y, p);
  465.     s = (fx - fy);
  466.     /* calculate the determinant */
  467.     if (s == 0.0)
  468.         r.x = MAXDOUBLE;
  469.     else
  470.         r.x = (-fx*(s - (fx*(x-y)*dpoly(y,p))))/(s*s);
  471.     /* calculate the trace later */
  472.     r.y = trace(x, y, p);
  473.     return(r);
  474. }
  475.  
  476. /* From the book "Symmetry In Chaos" by Michael Field and Martin Golubitsky.
  477.  * Written as a map of the real plane :
  478.  *     F(x,y) = m(x,y) + (u,v) + A(sin(2*Pi*x), sin(2*Pi*y))
  479.  *            + B(sin(2*Pi*x)*cos(2*Pi*y), sin(2*Pi*y)*cos(2*Pi*x))
  480.  *            + C(sin(4*Pi*x), sin(4*Pi*y))
  481.  *            + D(sin(6*Pi*x)*cos(4*Pi*y), sin(6*Pi*y)*cos(4*Pi*x))
  482.  *            + E(sin(2*Pi*y), sin(2*Pi*x))
  483.  * In addition, i allow the use of a different set of parameters for x and y.
  484.  * where A, B, C, D, E, m, u and v are all real parameters.
  485.  */
  486. pair
  487. golub(x, y, p)
  488. double x, y;
  489. double *p;
  490. {
  491.     pair z;
  492.     static double sx, sy, cy, cx, fpx, fpy;
  493.     extern double sin(), cos();
  494.  
  495.     sx = sin(M_2PI * x); sy = sin(M_2PI * y);
  496.     cx = cos(M_2PI * x); cy = cos(M_2PI * y);
  497.     fpx = 2.0 * M_2PI * x; fpy = 2.0 * M_2PI * y;
  498.     z.x = (p[6]*x) + p[5] + (p[0]*sx) + (p[1]*sx*cy) + (p[2]*sin(fpx));
  499.     z.x += ((p[3]*sin(3.0*M_2PI*x)*cos(fpy)) - (p[4]*sy));
  500.     z.y = (p[18]*y) + p[17] + (p[12]*sy) + (p[13]*sy*cx) + (p[14]*sin(fpy));
  501.     z.y += ((p[15]*sin(3.0*M_2PI*y)*cos(fpx)) - (p[16]*sx));
  502.     return(z);
  503. }
  504.  
  505. pair
  506. dgolub(x, y, p)    /* Not Done Yet - use the "-numeric" argument */
  507. double x, y;
  508. double *p;
  509. {
  510.     extern pair dnumeric();
  511.  
  512.     return(dnumeric(x, y, p));
  513. }
  514.  
  515. /* From the book "Symmetry In Chaos" by Michael Field and Martin Golubitsky.
  516.  * With additional parameters "n" and "p", the exponent and power, as doubles.
  517.  * Written as a map of the complex plane :
  518.  *     F(z) = (A+B*z*zbar+C*Re(z^n)+F*Re((z/|z|)^np)*|z|+Ei)*z + D*(zbar^(n-1))
  519.  * where A, B, C, D, E, F, p and n are all real parameters, z =x+iy complex.
  520.  * This is a combination of their two "symmetric  icon" formulas. In this
  521.  * more general form, i've kept the omega term to allow Zn symmetry and added
  522.  * the non-polynomial term to allow introduction of a singularity at the origin.
  523.  * If you want to get the first symmetric icon formula from the book, then just
  524.  * set F=0 above. If you want the second one from the book, then set E=0.
  525.  */
  526. pair
  527. marty(x, y, p)
  528. double x, y;
  529. double *p;
  530. {
  531.     pair t, z, w;
  532.     static double s, c;
  533.     extern void cmul(), zbar(), zpow(), cadd();
  534.     extern double rpzn(), sqrt();
  535.  
  536.     if (p[6] == 0.0)
  537.         s = 0.0;
  538.     else {
  539.         c = sqrt(x*x+y*y);
  540.         if (c == 0.0)
  541.             s = 0.0;
  542.         else
  543.             s = p[6]*rpzn(x/c,y/c,p[5]*p[7])*c;
  544.     }
  545.     cmul(&z,p[0]+(p[1]*((x*x)+(y*y)))+(p[2]*rpzn(x,y,p[5]))+s,p[4],x,y);
  546.     zbar(&w, x, y);
  547.     zpow(&t, w.x, w.y, p[5] - 1.0);
  548.     t.x *= p[3]; t.y *= p[3];
  549.     cadd(&w, z.x, z.y, t.x, t.y);
  550.     return(w);
  551. }
  552.  
  553. pair
  554. dmarty(x, y, p)    /* Not Done Yet - use the "-numeric" argument */
  555. double x, y;
  556. double *p;
  557. {
  558.     extern pair dnumeric();
  559.  
  560.     return(dnumeric(x, y, p));
  561. }
  562.  
  563. pair
  564. rotor(x, y, p) /* the rotor map (see Yorke's dynamics software) */
  565. double x, y; 
  566. double *p;
  567. {
  568.     static pair coord;
  569.     static double s;
  570.     extern double sin();
  571.  
  572.     s = sin(wrap(x, 0.0, M_2PI));
  573.     if (p[1] == 0.0) {
  574.         coord.y = y;
  575.         coord.x = x;
  576.         return(coord);
  577.     }
  578.     coord.y = (y - (p[0] * s))/p[1];
  579.     coord.x = wrap((x - coord.y), -M_PI, M_PI);
  580.     if (p[1] == 1.0)
  581.         coord.y = wrap(coord.y, -M_PI, M_PI);
  582.     return(coord);
  583. }
  584.  
  585. pair
  586. drotor(x, y, p)
  587. double x, y;
  588. double *p;
  589. {
  590.     pair r;
  591.  
  592.     r.x = r.y = MAXDOUBLE;
  593.     /* calculate the determinant */
  594.     if (p[1] == 0.0)
  595.         return(r);
  596.     r.x = 1.0 / p[1];
  597.     r.y = trace(x, y, p);
  598.     return(r);
  599. }
  600.  
  601. pair
  602. twistandflip(x, y, p) /* the Twist and Flip map (see Brown and Chua) */
  603. double x, y; 
  604. double *p;
  605. {
  606.     static pair coord;
  607.     static double s, c, theta, xmina;
  608.     extern double cos(), sqrt(), sin();
  609.  
  610.     xmina = x - p[1];
  611.     theta = (M_PI / p[0]) * sqrt((xmina*xmina) + (y*y));
  612.     s = sin(wrap(theta, 0.0, M_2PI));
  613.     c = cos(wrap(theta, 0.0, M_2PI));
  614.     coord.x = -1.0 * ((xmina * c) - (y * s) + p[1]);
  615.     coord.y = -1.0 * ((xmina * s) + (y * c));
  616.     return(coord);
  617. }
  618.  
  619. pair
  620. dtwistandflip(x, y, p)
  621. double x, y;
  622. double *p;
  623. {
  624.     pair r;
  625.     static double s, c, theta, t, xmina, poverw, yoverr; 
  626.     static double dxdx, dxdy, dydy, dydx;
  627.     extern double cos(), sqrt(), sin();
  628.  
  629.     xmina = x - p[1];
  630.     poverw = M_PI / p[0];
  631.     t = sqrt((xmina*xmina) + (y*y));
  632.     yoverr = y / t;
  633.     theta = poverw * t;
  634.     s = sin(wrap(theta, 0.0, M_2PI));
  635.     c = cos(wrap(theta, 0.0, M_2PI));
  636.     dxdx = ((-1.0) * poverw * ((xmina / t) * ((y*c) - (xmina*s)))) - c;
  637.     dydx = ((-1.0) * poverw * ((xmina / t) * ((xmina*c) - (y*s)))) - s;
  638.     dxdy = (poverw * yoverr *(xmina*s + (y*c))) + s;
  639.     dydy = ((-1.0) * poverw * yoverr *(xmina*c + (y*s))) - c;
  640.     /* calculate the determinant and trace of the Jacobian */
  641.     r.x = (dxdx*dydy) - (dydx*dxdy);
  642.     r.y = dxdx + dydy;
  643.     return(r);
  644. }
  645.  
  646. double
  647. Fmira(x, p)
  648. double x, p;
  649. {
  650.     static double z;
  651.  
  652.     z = x*x;
  653.     return(p*x + ((2.0*(1.0-p)*z)/(1.0 + z)));
  654. }
  655.  
  656. pair
  657. mira(x, y, p)        /* From Transformations Ponctuelles et Leurs */
  658. double x, y;        /* Applications (1973) by Bernussou, Hsu & Mira */
  659. double *p;
  660. {
  661.     pair coord;
  662.  
  663.     coord.x= y + Fmira(x, p[0]) + (p[1]*y*(1 + (p[2]*y*y)));
  664.     coord.y= Fmira(coord.x, p[0]) - x;
  665.     return(coord);
  666. }
  667.  
  668. double
  669. Fdmira(x, p)
  670. double x, p;
  671. {
  672.     static double z;
  673.  
  674.     return((4.0*(1.0-p)*x)/((1.0+z)*(1.0+z)));
  675. }
  676.  
  677. pair
  678. dmira(x, y, p)
  679. double x, y;
  680. double *p;
  681. {
  682.     static double a, b, c, d;
  683.     pair coord, r;
  684.  
  685.     coord = mira(x, y, p);
  686.     a = p[0] + Fdmira(x, p[0]);
  687.     b = 1.0 + p[1] + (3.0*p[1]*p[2]*y*y);
  688.     c = (Fdmira(coord.x, p[0]) * Fdmira(x, p[0])) - 1.0;
  689.     d = p[0] + Fdmira(coord.x, p[0]);
  690.     /* calculate the determinant of the jacobian */
  691.     r.x = b*((a*d) - c);
  692.     /* calculate the trace later */
  693.     r.y = trace(x, y, p);
  694.     return(r);
  695. }
  696.  
  697. pair
  698. lorenz(x, y, p)        /* From Physica D 35 (1989) 299-317 by Lorenz */
  699. double x, y; 
  700. double *p;
  701. {
  702.     pair coord;
  703.  
  704.     coord.x= ((1.0 + (p[1]*p[0]))*x) - (p[0]*x*y);
  705.     coord.y= ((1.0 - p[0]) * y) + (p[0]*x*x);
  706.     return(coord);
  707. }
  708.  
  709. pair
  710. dlorenz(x, y, p)
  711. double x, y;
  712. double *p;
  713. {
  714.     static double a, b, c;
  715.     pair r;
  716.  
  717.     a = 1.0 + (p[0]*p[1]) - (p[0]*y);
  718.     b = p[0]*x;
  719.     c = 1.0 - p[0];
  720.     /* calculate the determinant */
  721.     r.x = (a*c) - (2.0*b*b);
  722.     /* calculate the trace */
  723.     r.y = a + c;
  724.     return(r);
  725. }
  726.  
  727. pair
  728. logistic(x, y, p)        /* the "double logistic" map */
  729. double x, y; 
  730. double *p;
  731. {
  732.     pair coord;
  733.  
  734.     coord.x=(4.0*p[0]*x*(1.0 - x)) + ((1.0 - p[0]) * y);
  735.     coord.y=(4.0*p[1]*y*(1.0 - y)) + ((1.0 - p[1]) * x);
  736.     return(coord);
  737. }
  738.  
  739. pair
  740. dlogistic(x, y, p)
  741. double x, y;
  742. double *p;
  743. {
  744.     pair r;
  745.     double A, B;
  746.  
  747.     A = p[0]; B = p[1];
  748.     /* calculate the determinant */
  749.     r.x = ((16.0*A*B)*(1.0 - (2.0*y))*(1.0 - (2.0*x))) - ((1.0-A)*(1.0-B));
  750.     /* calculate the trace */
  751.     r.y = (4.0*(p[0] + p[1])) - (8.0*((p[0]*x)+(p[1]*y)));
  752.     return(r);
  753. }
  754.  
  755. pair
  756. dorband(x, y, p)        /* the Dorband "double logistic" map */
  757. double x, y; 
  758. double *p;
  759. {
  760.     pair coord;
  761.  
  762.     coord.x=(4.0*p[0]*y*(1.0 - y)) + ((1.0 - p[0]) * x);
  763.     coord.y=(4.0*p[1]*x*(1.0 - x)) + ((1.0 - p[1]) * y);
  764.     return(coord);
  765. }
  766.  
  767. pair
  768. ddorband(x, y, p)
  769. double x, y;
  770. double *p;
  771. {
  772.     pair r;
  773.  
  774.     /* calculate the determinant */
  775.     r.x=((1.0-p[0])*(1.0-p[1]))-((16.0*p[0]*p[1])*(1.0-(2.0*y))*(1.0-(2.0*x)));
  776.     /* calculate the trace */
  777.     r.y = 2.0 - p[0] - p[1];
  778.     return(r);
  779. }
  780.  
  781. double
  782. fpei(x,y, p)
  783. double x, y; 
  784. double *p;
  785. {
  786.     return((p[3]*x) - (p[4]*x*y));
  787. }
  788.  
  789. double
  790. gpei(x, y, p)
  791. double x, y;
  792. double *p;
  793. {
  794.     return((p[6]*x*y) - (p[5]*y));
  795. }
  796.  
  797. double
  798. ppei(x,y,p)
  799. double x, y;
  800. double *p;
  801. {
  802.     return(x + (p[1]*fpei(x,y,p)));
  803. }
  804.  
  805. double
  806. qpei(x,y,p)
  807. double x, y;
  808. double *p;
  809. {
  810.     return(y + (p[2]*gpei(x,y,p)));
  811. }
  812.  
  813. pair
  814. peitgen(x, y, p)    /* From "The Beauty of Fractals" */
  815. double x, y;        /* by Peitgen et al , page 125 */
  816. double *p;
  817. {
  818.     pair coord;
  819.  
  820.     coord.x=x+((p[0]*0.5)*(fpei(x,y,p)+fpei(ppei(x,y,p),qpei(x,y,p),p)));
  821.     coord.y=y+((p[0]*0.5)*(gpei(x,y,p)+gpei(ppei(x,y,p),qpei(x,y,p),p)));
  822.     return(coord);
  823. }
  824.  
  825. pair
  826. dpeitgen(x, y, p)
  827. double x, y;
  828. double *p;
  829. {
  830.     static double hover2, pfx, pfy, pgx, pgy, a, b, c, d;
  831.     pair r;
  832.  
  833. printf("");
  834.     hover2 = p[0] * 0.5;
  835.     a = (p[2]*p[4]*p[5])-(2.0*p[1]*p[3]*p[4])+(p[1]*p[2]*p[3]*p[4]*p[5]);
  836.     b = (p[2]*p[4]*p[6]) + (p[1]*p[2]*p[3]*p[4]*p[6]);
  837.     c = p[1]*p[2]*p[4]*p[4]*p[5]; 
  838.     d = p[1]*p[2]*p[4]*p[4]*p[6];
  839.     pfx = p[3] - (p[4]*y) + (p[3]+(p[1]*p[3]*p[3])) + (a*y) - 
  840.             (2.0*b*x*y) - (c*y*y) + (2.0*d*x*y*y);
  841.     pfy = (p[1]*p[4]*p[4]) - (p[4]*x) + (a*x) - (b*x*x) - (2.0*c*x*y) + 
  842.             (2.0*d*x*x*y);
  843.     a=p[6]-(2.0*p[2]*p[6]*p[5])+(p[1]*p[3]*p[6])-(p[1]*p[2]*p[3]*p[6]*p[5]);
  844.     b = (p[2]*p[6]*p[6]) + (p[1]*p[2]*p[3]*p[6]*p[6]);
  845.     c = p[1]*p[2]*p[4]*p[6]*p[5]; 
  846.     d = p[1]*p[2]*p[4]*p[6]*p[6];
  847.     pgx = (p[6]*y) + (a*y) + (2.0*b*x*y) + (c*y*y) - (2.0*d*x*y*y);
  848.     pgy = (p[6]*x) - (2.0*p[5]) + (p[2]*p[5]*p[5]) + (a*x) + (b*x*x) + 
  849.             (2.0*c*x*y) - (2.0*d*x*x*y);
  850.     a = 1.0 + (hover2 * pfx);
  851.     b = hover2 * pfy;
  852.     c = hover2 * pgx;
  853.     d = 1.0 + (hover2 * pgy);
  854.     /* calculate the determinant of the jacobian */
  855.     r.x = (a*d) - (b*c);
  856.     /* calculate the trace */
  857.     r.y = a + d;
  858.     return(r);
  859. }
  860.  
  861. pair
  862. julia(x, y, p)        /* the Julia-Mandelbrot map using real rather than */
  863. double x, y;         /* complex coordinates */
  864. double *p;
  865. {
  866.     pair coord;
  867.  
  868.     coord.x = (x*x) - (y*y) + p[0] + (p[2]*x);
  869.     coord.y = (2.0*x*y) + p[1] + (p[3]*p[2]*y);
  870.     return(coord);
  871. }
  872.  
  873. pair
  874. djulia(x, y, p)
  875. double x, y;
  876. double *p;
  877. {
  878.     pair r;
  879.  
  880.     /* calculate the determinant */
  881.     r.x = (4.0*((x*x) + (y*y))) + (p[2]*((2.0*x*(p[3]+1)) + (p[2]*p[3])));
  882.     /* calculate the trace */
  883.     r.y = (4.0*x) + p[2] + (p[2]*p[3]);
  884.     return(r);
  885. }
  886.  
  887. pair
  888. brussel(x, y, p)    /* The Brusselator flow after numerically integrating */
  889. double x, y;        /* using the Euler method with time step p[0] */
  890. double *p;
  891. {
  892.     pair coord;
  893.     static double a, b;
  894.  
  895.     a = x*x*y;
  896.     b = p[1]*x;
  897.     coord.x = x + (p[0] * (p[2] - b + a - x));
  898.     coord.y = y + (p[0] * (b - a));
  899.     return(coord);
  900. }
  901.  
  902. pair
  903. dbrussel(x, y, p)
  904. double x, y;
  905. double *p;
  906. {
  907.     double a, b, c, d;
  908.     pair r;
  909.  
  910.     d = 2.0*p[0]*x*y;
  911.     a = 1.0 - (p[0]*p[1]) + d - p[0];
  912.     b = p[0]*x*x;
  913.     c = (p[0]*p[1]) - d;
  914.     /* calculate the determinant */
  915.     d = 1.0 - b;
  916.     r.x = (a*d) - (b*c);
  917.     /* calculate the trace */
  918.     r.y = a + d;
  919.     return(r);
  920. }
  921.  
  922. pair
  923. baru(x, y, p)    /* the Barugola map   "sur certaines zones chaotiques*/
  924. double x, y; 
  925. double *p;
  926. {
  927.     pair coord;
  928.     static double A;
  929.     extern double exp();
  930.  
  931.     A = ((p[3] * p[0]) / (1.0 - exp(p[2] * p[0])));
  932.     coord.x = (x * exp(p[0] * (1.0 - p[1] * x) - (A * y)));
  933.     coord.y = (x * ( 1.0 - exp(-A * y)));
  934.     return(coord);
  935. }
  936.  
  937. pair
  938. dbaru(x, y, p)
  939. double x, y;
  940. double *p;
  941. {
  942.     static pair r;
  943.     static double A;
  944.     extern double exp();
  945.  
  946.     A = ((p[3] * p[0]) / (1.0 - exp(p[2] * p[0])));
  947.     /* calculate the determinant */
  948.     r.x = A*x*(1.0 - p[1] * p[0] * x) * exp((p[0] * (1.0 - p[1] * x)-(A * y)));
  949.     /* calculate the trace later */
  950.     r.y = trace(x, y, p);
  951.     return(r);
  952. }
  953.  
  954. #ifdef Plendo
  955. pair
  956. q1(x, y, p)    /* the Mira Quad map  */
  957. double x, y; 
  958. double *p;
  959. {
  960.     pair coord;
  961.     static double xsq, ysq;
  962.  
  963.     xsq = x * x; ysq = y * y;
  964.     coord.x = (p[0]*x) + (p[2]*y) + p[3] + (p[4]*xsq) + (p[5]*ysq) + p[6];
  965.     coord.y = (p[7]*x) + (p[8]*y) + p[1] + (p[9]*xsq) + (p[10]*ysq) + p[11];
  966.     return(coord);
  967. }
  968.  
  969. pair
  970. dq1(x, y, p)
  971. double x, y;
  972. double *p;
  973. {
  974.     pair r; 
  975.     double pfx, pgy;
  976.  
  977.     pfx = p[0] + (2.0*p[4]*x);
  978.     pgy = p[8]+(2.0*p[10]*y);
  979.     r.x = pfx * pgy - ((p[7] + (2.0*p[9]*x))*(p[2]+(2.0*p[5]*y)));
  980.     r.y = pfx + pgy;
  981.     return(r);
  982. }
  983.  
  984. pair
  985. q2(x, y, p)    /* the Mira Agg map  */
  986. double x, y; 
  987. double *p;
  988. {
  989.     pair coord;
  990.  
  991.     coord.x = (p[0]*x) + (p[2]*y);
  992.     coord.y = (p[3]*x) + (p[5]*y) + + (p[4]*x*x) + (p[1]*x*y);
  993.     return(coord);
  994. }
  995.  
  996. pair
  997. dq2(x, y, p)
  998. double x, y;
  999. double *p;
  1000. {
  1001.     pair r;
  1002.  
  1003.     r.x = p[0]*(p[5]+(p[1]*x));
  1004.     r.x -= p[2]*(p[3]+(p[1]*y)+(2.0*p[4]*x));
  1005.     r.y = p[0] + p[5] + (p[1]*x);
  1006.     return(r);
  1007. }
  1008.  
  1009. pair
  1010. q3(x, y, p)    /* the Mira-Narayaninsamy map  */
  1011. double x, y; 
  1012. double *p;
  1013. {
  1014.     pair coord;
  1015.  
  1016.     coord.x = p[0] + (p[1]*x) + (x*x) - (y*y);
  1017.     coord.y = (2.0*x*y) - (2.5*p[1]*y);
  1018.     return(coord);
  1019. }
  1020.  
  1021. pair
  1022. dq3(x, y, p)
  1023. double x, y;
  1024. double *p;
  1025. {
  1026.     pair r;
  1027.  
  1028.     r.x = 4.0*((x*x) + (y*y));
  1029.     r.x -= 3.0*p[1]*x;
  1030.     r.x -= 2.5*p[1]*p[1];
  1031.     r.y = (4.0*x) - (1.5*p[1]);
  1032.     return(r);
  1033. }
  1034.  
  1035. pair
  1036. c1(x, y, p)    /* the Kawakami-Kawasaki map  */
  1037. double x, y; 
  1038. double *p;
  1039. {
  1040.     pair coord;
  1041.  
  1042.     coord.x = (p[1]*x) + (p[2]*y);
  1043.     coord.y = (p[0]*x) + (p[3]*y) + (p[4]*x*x*x);
  1044.     return(coord);
  1045. }
  1046.  
  1047. pair
  1048. dc1(x, y, p)
  1049. double x, y;
  1050. double *p;
  1051. {
  1052.     pair r;
  1053.  
  1054.     r.x = p[1]*p[3];
  1055.     r.x -= p[2]*(p[0] + (3.0*p[4]*x*x));
  1056.     r.y = p[1] + p[3];
  1057.     return(r);
  1058. }
  1059. #endif Plendo
  1060.  
  1061. pair
  1062. dnumeric(x, y, p)
  1063. double x, y;
  1064. double *p;
  1065. {
  1066.     static pair c, cx, cy;
  1067.     static pair r;
  1068.     extern PFP map;
  1069.     extern double delta;
  1070.  
  1071.     c = (*map)(x, y, p);
  1072.     cx = (*map)(x+delta, y, p);
  1073.     cy = (*map)(x, y+delta, p);
  1074.  
  1075.     /* first, calculate the determinant of the Jacobian */
  1076.     r.x = (cx.x*cy.y)-(cx.x*c.y)-(c.x*cy.y)-(cy.x*cx.y)+(cy.x*c.y)+(c.x*cx.y);
  1077.     r.x /= (delta*delta);
  1078.     /* next, calculate its trace */
  1079.     r.y = cx.x + cy.y - c.x - c.y;
  1080.     r.y /= delta;
  1081.     return(r);
  1082. }
  1083.